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Abstract 

Travel-time  tomography  was  vised  to  obtain  the  lateral  variations  in  Pn  veloci¬ 
ties  and  Pn  velocity  azimuthal  anisotropy  beneath  the  Basin  and  Range  Province.  We 
used  1226  Pn  arrivals  from  184  earthquakes  recorded  by  113  stations.  The  tomographic 
image  shows  the  lateral  variation  in  Pn  velocity  from  7.9  km^s  in  the  northern  Province 
to  7.6  kmyfe  in  the  southern  Province.  Azimuthal  variations  in  Pn  velocities  are  less 
than  one  percent  as  determined  from  the  travel  time  residuals  after  taking  out  the  con¬ 
tribution  of  the  lateral  velocity  variations  of  the  Pn  tomography.  Although  the  standard 
errors  of  the  azimuthal  variations  are  about  two  percent,  there  is  no  reason  to  believe 
the  existence  of  the  Pn  azimuthal  anisotropy  in  the  uppermost  mantle  beneath  the  Basin 
and  Range  Province.  Beghoul  and  Barazangi’s  (1990)  suggestion  of  more  than  three 
percent  azimuthal  anisotropy  is  the  result  of  the  lateral  heterogeneity  and  mantle  velo¬ 
city  gradient  Before  tomographic  inversion,  we  refined  the  Pn  travel  times  by  consid¬ 
ering  the  lateral  variations  in  the  crustal  thickness  and  velocities  in  the  region  based  on 
the  numerical  experiment  and  theoretical  development  given  in  this  paper.  These 
refinements  are  crucial  in  obtaining  an  accurate  velocity  image.  The  average  velocity  is 
7.72  ±0.16  kmyfeec,  and  an  average  vertical  mantle  velocity  gradient  is  8.0  X10-1  s-1. 
The  resolution  for  Pn  velocity  image  is  generally  within  two  to  five  degrees. 


Introduction 

Anisotropy,  a  characteristics  of  many  mantle  minerals,  is  important  to  seismolo¬ 
gy,  since  it  is  responsible  for  the  large  variations  in  seismic  velocities.  Some  research¬ 
ers  have  suggested  seismic  velocity  variations  due  to  anisotropy  are  potentially  much 
greater  than  those  due  to  other  effects  such  as  composition  and  temperature.  Thus,  it  is 


necessary  to  understand  seismic  anisotropy  well  before  one  attempts  to  infer  chemical, 
mineralogical  and  temperature  variations  from  seismic  data  [e.  g.,  Anderson,  1889]. 
However,  some  evidences  of  the  existence  of  seismic  anisotropy  can  be  interpreted  as 
the  lateral  heterogeneity  of  the  structure,  since  there  is  often  a  direct  trade-off  between 
anisotropy  and  heterogeneity.  In  order  to  circumvent  this  problem,  we  present  a 
method  of  identifying  Pn  velocity  azimuthal  anisotropy  in  the  presence  of  the  lateral 
heterogeneity  of  the  crustal  and  upper  mantle. 

An  anisotropic  upper  mantle  would  be  expected  if  it  is  olivine  rich  [  e.  g., 
Duffy  and  Anderson,  1989;  p.  306,  Anderson,  1989  ].  Crystal  alignment  may  be 
caused  by  the  motions  in  the  mantle  associated  with  mantle  convection,  internal  defor¬ 
mation  of  the  lithosphere,  the  preferred  local  stress  field  orientation,  or  recrystalliza- 
tion.  Thus,  attempts  at  finding  evidence  of  seismic  anisotropy  are  largely  focused  on 
tectonically  active  regions  in  the  world.  Due  to  limited  data  and  complexity  of  the 
Earth,  especially  in  tectonically  active  regions,  it  is  difficult  to  identify  anisotropy 
unambiguously  in  the  presence  of  structural  heterogeneity.  In  order  to  minimize  the 
contribution  of  heterogeneity,  researchers  have  turned  to  shear  wave  splitting  to  toy  to 
find  anisotropy,  [  t.  g.,  Savage  et  al.  1990;  Silver  and  Chan,  1988,  1991;  Xie,  1992  ]. 
Numerical  experiments  show  that  coupled  mode  waveform  anomalies  of  surface  waves 
in  the  frequency  range  of  0-20  mHz  can  be  diagnostic  of  upper  mantle  azimuthal  an¬ 
isotropy  [  Park  and  Yu,  1992;  Tanimoto  and  Anderson,  1985;  Yu  and  Park,  1992  ]. 

Pn  travel  times  have  been  used  to  find  evidence  for  the  existence  of  azimuthal 
anisotropy  [  e.  g.,  Bamford  et  al.,  1979;  Begboul  and  Barazangi,  1990;  Hess,  1964; 
Rait  et  al.,  1971;  Vetter  and  Minster,  1981  ].  For  the  Basin  and  Range  province, 
Beghoul  and  Barazangi  [  1990  ]  suggested  a  3.2  percent  azimuthal  anisotropy  in  Pn 
velocity.  In  this  study,  we  rigorously  examine  Pn  velocities  beneath  the  Basin  and 
Range  province.  We  used  all  possible  ISC  Pn  travel  times  from  1964  to  1990.  We  ac¬ 
count  for  the  three  dimensional  local  crustal  and  mantle  velocity  structure  (  including 
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the  topography  of  the  surface  and  the  Moho  )  to  refine  the  apparent  Pn  velocities  be¬ 
fore  inversion  employing  the  numerical  method  discussed  in  detail  later.  After  tomo¬ 
graphic  inversion,  we  allocated  the  non- diminished  Pn  travel  time  residuals  back  to  the 
vertical  travel  times  of  the  source  and  receiver.  Ibis  process  is  actually  the  source  re¬ 
location  and  local  velocity  structure  refinement  After  this  process,  we  obtain  a  new 
data  set  of  apparent  Pn  velocities.  We  then  invert  the  new  data  set  and  allocate  the 
non-diminished  residuals  after  inversion  again  until  the  allocation  so  small  that  the 
results  of  the  last  two  iterations  are  almost  the  same.  In  this  study  we  used  four  itera¬ 
tions.  After  the  whole  process  is  done,  we  assume  that  the  non-diminished  Pn  travel 
times  residuals  are  mainly  due  to  the  azimuthal  velocity  dependence,  since  the  contri¬ 
bution  of  lateral  variation  of  the  Pn  velocities  are  taken  out  by  the  tomographic  inver¬ 
sion.  Thus,  we  can  find  out  if  Pn  velocity  residuals  calculated  from  travel  time  residu¬ 
als  depend  on  azimuth. 


Data 

Using  the  ISC  database,  we  collected  1226  Pn  travel  times  from  184  earth¬ 
quakes  (  Ml  >2.9  )  that  occurred  in  or  around  the  Basin  and  Range  province  from 
1964  to  1990,  recorded  at  113  stations.  The  ray  paths  are  shown  in  Figure  1.  We 
selected  paths  for  which  the  apparent  Pn  velocity  is  between  7.0  -  8.5  kmyfeec,  after 
correcting  the  average  crustal  thickness  of  35  km  [  Braile  et  a/.,  1989  ].  The  distance 
range  is  2-11  degrees. 


Method  and  Procedure 

We  begin  with  the  observed  Pn  travel  times,  t^,  which  can  be  written  as: 


W  -or+0+7+JSdl,  (1) 

S 


where  a  is  the  vertical  travel  time  from  source  to  Moho,  /?  is  that  from  Moho  to 
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receiver,  and  'y  is  the  contribution  of  the  vertical  velocity  gradient  of  the  mantle.  The 

r 

horizontal  travel  time,  J Sdl  is  from  source,  s,  to  receiver,  r,  where  S  is  local  slowness. 

8 

In  the  following,  we  give  a  detailed  description  how  to  estimate  a,  /?,  and 

Two  dimensional  crustal  velocity  structure 

In  a  crustal  structure  with  weak  heterogeneity,  a  and  can  be  written  approxi¬ 
mately  as: 

j-c,  i-c, 

and  (2) 

H  i-l 

where  and  iij.  are  the  numbers  of  layers  beneath  the  source  and  the  receiver  to  the 
local  Moho,  respectively.  Hj  is  the  layer  thickness  beneath  the  source,  while  H,  is  the 
layer  thickness  beneath  the  receiver,  r]l  “(l^,2  —  p2)1/2  is  vertical  slowness,  where  Vj 
is  velocity,  and  p  is  ray  parameter,  a  is  calculated  from  the  velocity  structure  close  to 
the  source,  and  /?  is  from  the  structure  close  to  the  receiver  (  the  topography  of  the 
surface  can  be  easily  included  ). 

We  can  show'  the  validity  of  Eq.  (2)  by  comparing  with  finite  difference  syn¬ 
thetic  seismograms  (  Vidale  et  ah,  1985  ).  Figure  2  gives  a  two  dimensional  model, 
for  which  the  crustal  thickness  varies  from  50  km  to  30  km.  Figure  3  gives  the  finite 
difference  synthetics  of  the  model  (  solid  lines,  which  are  truncated  seismograms  ), 
predictions  of  Eq.  (2)  (  arrows  )  ,  and  picks  (  with  0.1  seconds  error  bar  )  from  the 
synthetics  (  solid  dots  ).  The  maximum  difference  between  picked  and  calculated 
travel  times  is  0.3  seconds,  which  contrast  greatly  with  the  more  than  one  second 
difference  in  the  reduced  travel  times  for  all  traces.  Assuming  the  model  at  the  source 
end  is  the  model  for  the  whole  path,  we  have  a  reduced  travel  time  of  7.9  seconds.  In 
contrast,  the  reduced  travel  time  is  only  4.6  seconds  if  using  the  model  at  200  km  from 
tiie  source  as  the  model  for  the  whole  path.  Thus,  Eq.  (1)  is  a  good  approximation  for 
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the  head  travel  times  in  two  dimensional  velocity  structures.  If  the  crustal  velocity 
structure  does  not  change  very  rapidly  in  a  few  tens  of  kilometers  (  order  of  the  crustal 
thickness  )  in  horizontal  direction,  the  error  of  the  predictions  of  Eq.  (2)  is  very  small  ( 
see  the  210-270  km  and  360-470  km  ranges  of  Figures  2  and  3  ). 

The  model  used  to  estimate  a  is  the  model  at  50  km  from  the  source,  roughly 
the  distance  from  the  source  to  Moho.  For  /#,  the  model  used  is  at  40  km  to  the 
receiver.  If  the  velocity  structure  is  roughly  the  same  in  a  few  tens  of  kilometers  in  the 
horizontal  direction,  which  is  not  true  for  the  model  given  in  Figure  2,  we  can  use  the 
velocity  structure  below  the  source  and  receiver  to  estimate  a  and  /?  of  Eq.  (2). 

The  absolute  travel  times  of  finite  difference  synthetic  seismograms  are  refined 
by  comparing  the  finite  difference  results  of  an  averaged  one- dimensional  model  of  the 
a  two-dimensional  model  used  (Figure  2).  For  one-dimensional  models,  the  head  wave 
travel  times  can  be  calculated  analytically.  The  highest  frequency  end  of  the  synthetics 
is  one  Hertz.  In  the  complex  structures,  the  picked  P„  time  may  not  be  the  minimum 
travel  time  for  the  model,  because  the  energy  of  the  phase  may  be  too  small  to  be 
noticed. 

The  good  approximation  of  Eqs.(l)  and  (2)  can  be  understood  theoretically.  The 
head  waves  travel  almost  horizontally  in  the  mantle  so  that  the  horizontal  travel  time 
would  be  approximately  the  integral  in  Eq.  (1).  The  head  wave  travel  time  would  be 
the  minimum  travel  time,  so  the  ray  has  the  critical  angle  incidence  to  the  local  Moho 
at  the  source  and  leaves  the  local  Moho  at  the  critical  angle  to  go  to  the  receiver. 
Eq.(l)  is  a  good  approximation  as  long  as  the  structure  'is  sufficiently  close  to  one¬ 
dimensional,  since  the  ray  parameter  P  is  no  longer  a  constant  in  two-dimensional 
velocity  structures. 

A  constant  velocity  gradient  mantle 

If  a  mantle  has  a  constant  velocity  gradient,  then  the  travel  times  of  head  waves 
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are  not  a  linear  function  of  distance.  If  we  use  the  form 

v-'Voli+c2)  (3) 

to  represent  the  velocities  in  the  mantle,  where  z  is  the  depth  from  Moho  (  downward 
is  positive  ),  c  is  a  constant  whose  product  with  v0  gives  the  mantle  velocity  gradient 
for  a  flat  Earth.  For  a  spherical  Earth,  the  velocity  gradient  is  (  c  —1.58  XIO"4)^ 
when  the  units  for  c  and  v0  are  km-1  and  kmyfe,  respectively  [Helmberger,  1973].  The 
observed  bead  wave  velocity  is  vob6  /i,  where  X  is  the  horizontal  distance  in  the 
mantle,  and  t  is  the  minimum  travel  time  of  the  model  at  this  distance  (  we  used  vobs 
here  because  its  definition  is  the  same  as  that  used  to  estimate  the  average  observed 
velocity  for  head  waves,  and  to  show  the  difference  from  the  conventionally  used  head 
wave  velocity  of  v0).  vote  can  be  written  as 

vob6~vo+-^-voX2,  (4) 

when  c  Cl,  which  is  usually  holds  for  the  real  Earth.  Eq.  (4)  is  valid  when 

X  <Xm  «((1  -f  cHm)2— 1)1/2  /c,  where  Hm  is  the  maximum  depth  that  Equation  (3) 

holds.  When  X  >Xm,  we  have 

2X 

vob6  **0  +vocHm(l — g“)  (5) 

The  derivations  of  Eqs.  (4)  and  (5)  are  given  in  the  Appendix. 

From  Eqs.  (4)  and  (5),  7  in  Eq.  (1)  can  be  written  as 

_  X  6v 

7  - - , 

vo  vo 

where  6v  is  the  second  term  in  Eqs.  (4)  or  (5). 

Inversion  method 

Once  a  and  /?  are  calculated  from  Eq.  (2)  and  the  local  velocity  structure  is 
removed  from  in  Eq.  (1),  we  have 
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Voi* S3')  and  X -D -A -B,  (7) 

where  D  is  the  horizontal  source-receiver  distance,  the  primed  quantities  are  estimated 
ones  from  available  information.  A  and  B,  the  horizontal  travel  distances  of  the  PD  ray 
in  the  crust  beneath  the  source  and  receiver,  can  be  estimated  from  local  crustal  model. 
For  every  path,  once  vob6  and  X  are  determined,  we  can  calculate  c  from  Eqs.  (4)  and 
(5)  using  the  least  squared  method,  and  then  7  from  Eq.  (6). 

For  each  of  ( i ,  j ),  source  i  to  receiver  j,  we  have: 

(J S  dl)i,  -B^ijkSk  -(W)ij-°  W '-•»'«  -RES*  (8) 

8 

where  is  the  length  of  the  line  segment  over  which  the  (i  j)th  path  overlaps  the 
kth  grid,  sk  is  the  slowness  in  the  kth  grid,  and  the  travel  time  residual  for  (i,  j)  path  is 

RESijWiai+^+^j 

'd+tfrP  'jMm'a)  • 

Then  Eq.  (8)  can  be  rewritten  as 

|>ijksk+REs  W  'r  ^  V  (0) 

The  right  side  of  Eq.  (9)  is  known.  This  equation  forms  a  standard  linear  inverse  prob¬ 
lem  for  unknown  sk.  We  used  a  back-projection,  or  SIRT  algorithm  (  e.  g.  Humphreys 
and  Clayton,  1998;  Xie  and  Mitchell,  1990  )  to  solve  for  sk.  In  particular,  we  use  the 
9-grid  smoothing  operation  (  Suetsugu  and  Nakanishi,  1985;  Xie  and  Mitchell,  1990  ) 
between  iterations.  This  operation  is  used  to  overcome  insufficient  spatial  data  cover¬ 
age,  at  a  slight  cost  to  resolution. 

Method  to  refine  o  and 

In  the  back-projection  calculation,  sk  values  are  updated  during  the  iterations  to 
reduce  the  misfit  of  left  and  right  side  of  Eq.  (9).  When  the  misfit  stops  decreasing, 
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we  end  up  with  with  the  final  slowness  model  (  values  of  sk  )  and  some  non- 
diminished  RESjj  values.  To  the  first  order  approximation  we  assume  that  RES,j  are 
dominatedly  caused  by  imprecise  estimates  on  a,  and  f3};  we  then  have 

RES;j  «5a1+^J.  (10) 

Non-zero  6a  are  primarily  due  to  the  mislocadon  of  the  source  (  location,  depth,  and 
origin  time  ),  and  the  erroneous  crustal  structure  used  before  the  inversion,  and  non¬ 
zero  6p  is  purely  from  erroneous  crustal  structure.  Eq.  (10)  also  forms  a  linear  inver¬ 
sion  problem,  for  which  we  used  the  least  square  method  to  allocate  REStJ  back  to  ax 
and  This  means  that  the  sources  are  relocated  and  the  local  crustal  models  for  the 
source  and  receiver  regions  have  been  improved  in  the  new  slowness  model. 

After  new  a,  and  ^  are  determined,  we  implement  the  back- projection  again  to 
update  sk  model,  then  update  a  and  again,  until  the  maximum  of  the  allocations  6a 
and  6fi  of  a  and  /?  reaches  some  small  limit 


Method  of  estimating  azimuthal  anisotropy 

Since  the  velocity  image  (  sk  model  )  is  an  average  in  all  directions,  it  is  rea¬ 
sonable  to  assume  that  the  non-zero  travel  time  residuals,  RES,j,  result  mainly  from  the 
possible  azimuthal  anisotropy,  after  all  inversion  process  in  the  earlier  sections  is  done. 
We  define  the  velocity  residual  as: 


6\  vob6  vsyn  ^  RES 
v  vsyn  W  —a -£-7' 


For  one  particular  azimuth  6,  we  use  average 
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as  the  velocity  residuals  at  this  azimuth,  where  i  is  the  path  within  the  20  degree 
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interval.  And  we  use  the  square  root  of 
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as  the  standard  deviation  of  tins  velocity  residual. 


Results 

We  used  the  resulting  models  of  seismic  sounding  studies  [  Braile  et  ah,  1989  ] 
as  local  models  to  estimate  a  and  0,  and  obtained  a '  ,  0 '  values  for  the  source  and 
the  receiver  delays,  calculated  from  Eq.  (2  ).  In  this  estimation,  we  assumed  that  the 
velocity  structure  does  not  vary  rapidly.  The  topography  of  the  surface  of  the  Earth  is 
included. 


Mantle  velocity  gradient 

Using  Burdick  and  Helmberger’s  [  1978  ]  T7  model,  which  was  proposed  for 
the  western  United  States,  we  can  get  an  estimate  for  Xg*  of  2400  km,  much  greater 
than  the  maximum  distance  we  used  in  this  study.  Thus,  we  can  use  Eqs.  (4)  and  (7) 
to  estimate  c,  and  obtain  c  *2.56  X10-4,  which  is  the  average  for  all  ray  paths,  and 
mantle  velocity  gradient  2.0  XlO^s-1,  including  the  contribution  of  the  Earth’s  spheri¬ 
city  1.2  X10"*s~*.  7'  can  be  estimated  by  using  Eq.  (6).  At  a  distance  of  1200  km, 
7'  is  -0j6  seconds. 

Pn  velocity  image 

After  accounting  for  a  \  0'  and  7',  we  used  four  allocation  iterations  to  find 
solutions  for  Eqs.  (9)  and  (10)  so  that  the  maximum  of  6a  and  60  was  less  than  0.2 
seconds.  We  divide  the  area  between  29°  -47®N  and  107®-126®W  into  1406  blocks  of 
a  size  of  0.5®  X0.5®  .  Figure  4  gives  the  final  velocity  image.  The  average  velocity  of 
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this  image  is  7.72  ±0.16  km^.  (The  Pn  velocities  are  given  at  depth  of  30  km.)  The 
velocity  image  shows  a  faster  velocity  in  the  north  and  a  slower  velocity  in  the  south  ( 
Figure  4  ).  The  slow  velocity  regions  in  the  uppermost  mantle  beneath  Sierra  Madre 
and  the  southern  end  of  Nevada  appear  to  be  connected. 

Along  the  Sierra  Nevada,  the  Pn  velocity  is  less  than  7.7  kmyfeec,  which  agrees 
well  with  the  velocity  7.6-7.65  km^  of  Jones  et  at.  |  1092  ],  7.6-7.7  km/fe  of  Hearn  et 
al.  [  1901  ]. 

The  slow  Pn  velocities  in  the  southern  Basin  and  Range  province  (Figure  4)  are 
compatible  with  Biasi  and  Humphreys’  (  1992  ]  velocity  image  of  depth  range  30-60 
km.  However,  they  obtained  a  faster  region  at  the  northwest  of  Southwest  Nevada  Vol¬ 
canic  Field  (  roughly  at  38°N,  117°W)  closer  to  the  Nevada  boundary  than  this  study. 
The  reason  for  this  is  probably  that  the  mantle  velocity  gradient  is  larger  beneath  that 
region.  Hearn  et  al.  [1991]  also  found  a  slow  region  in  the  south,  albeit  a  much 
smaller  one,  which  may  be  because  of  the  flat  Moho  assumption  used  in  their  study. 
We  will  return  to  this  later. 

Azimuthal  anisotropy 

The  analysis  of  the  travel  time  residuals,  after  the  inversion  and  allocation  pio- 
cedure  is  done,  is  given  in  Figure  5,  in  five  degree  spacing,  using  Eqs.  (12)  and  (13). 

#-H0 

The  denominators,  £(6  -  \B{  -6  |/2),  in  Eqs.  (12)  and  (13)  are  given  in  Table  1. 

#-40 

The  value  of  the  denominator  at  one  azimuth  represents  the  number  of  Pn  travel  time 
picks  in  a  20  degree  azimuthal  interval  centered  at  the  azimuth.  The  average  is  241.1, 
roughly  equivalent  to  70  picks  in  20  degree  intervals.  None  of  the  average  velocity 
residuals  exceeds  one  percent  (  Figure  5).  The  average  velocity  residuals  at  azimuth 
80°-  100°  are  about  0.0  percent  However,  the  denominator  values  for  their  20  degree 
interval  are  less  than  half  of  the  average  (  Table  1  ).  The  error  bars  are  mostly  larger 
than  two  percent  Thus,  the  azimuthal  anisotropy  in  Pn  velocity,  if  any,  is  not  likely  to 
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be  greater  than  one  percent  for  the  Basin  and  Range  province. 

Error  analysis  and  resolution  of  the  PD  velocity  image 

In  this  section,  we  quantitatively  assess  the  quality  of  the  two-dimensional  Pn 
velocity  image  obtained  in  this  study.  We  also  present  the  Pn  velocity  images,  (1) 
without  correcting  for  a  mantle  velocity  gradient  and  (2)  assuming  a  homogeneous 
crust  and  mantle. 

In  tomography,  the  spatial  resolution  at  a  given  grid  can  be  approximated  by 
the  point  spreading  function  (  p.s.f.  )  calculated  for  the  grid,  t.  g.,  Humphreys  and 
Clayton,  1988;  Xie  and  Mitchell,  1990.  We  calculated  the  p.s.f.  for  six  grids,  five  close 
to  the  boundaries,  one  in  the  center  of  the  region  (  Figure  6a  ).  We  put  three  grids  of 
the  p.s.f.  on  one  plot,  since  they  do  not  overlap  each  other.  The  spatial  resolution  is 
generally  less  than  3  degrees,  and  less  than  2  degrees  at  locations  where  the  ray  cover¬ 
age  is  more  dense.  However  the  grid  beneath  Sierra  Nevada  (  the  lower  dark  one  of 
Figure  6a  )  spreads  about  five  degrees.  This  is  caused  by  parallel  rays  passing  through 
most  of  area.  However,  point  spreading  function  gives  only  the  resolution  that  can  be 
achieved  from  the  given  ray  path  coverage. 

Figure  7  gives  the  velocity  image  without  correcting  for  mantle  velocity  gra¬ 
dient,  that  is,  7'  *0.0.  The  average  velocity  is  7.73  ±0.16  kmyfeec.  The  main  pattern 
of  the  velocity  distribution  of  Figure  7  is  about  same  as  that  of  Figure  4.  However,  the 
slower  Pn  velocity  regions  in  the  south  are  smaller  and  not  connected  to  each  other  in 
Figure  7,  as  compared  with  that  in  Figure  4.  The  small  difference  in  Figures  4  and  7  is 
due  to  small  positive  mantle  velocity  gradient,  2.0  XlO^s-1. 

As  has  been  discussed  before,  we  used  the  crustal  and  upper  mantle  velocity 
models  of  seismic  sounding  studies  to  estimate  a  '  and  /? '  from  Eq.  (2).  In  this  sec¬ 
tion,  we  show  the  results  from  the  assumption  of  a  homogeneous  uppermost  mantle 
velocity  model,  which  is  an  average  model  of  the  Basin  and  Range  province  [  Braile  et 
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al.,  1089  j,  with  a  compressions]  velocity  of  6.25  kmyfe  for  the  crust  and  7.8  kmyfe  for 
the  mantle,  and  with  a  thickness  of  35  km,  for  the  whole  Basin  and  Range  province. 
The  a '  and  /?'  values  estimated  also  from  Eq.  (2).  The  velocity  image  is  given  in  Fig¬ 
ure  8.  The  average  Pn  velocity  of  7.77  ±0.16  kmyfe  in  Figure  8  is  higher  than  that  of 
7.72  kmyfe  in  Figure  4.  The  main  features  of  the  velocity  distribution  are  similar 
except  for  the  southern  end,  at  which  the  Pn  velocities  are  much  exaggerated.  The  slow 
regions  are  reduced  beneath  the  Sierra  Nevada  and  southern  Nevada  (  Figures  4  and  8 
).  The  fast  regions  are  expanded,  for  example,  see  northern  portion  of  the  study  area. 
The  southern  end  of  the  study  area,  close  to  the  Gulf  of  California,  Pn  velocities 
becomes  fast  (  Figure  8  ),  which  is  veiy  similar  to  the  results  of  Hearn  et  al.  [  1991  ]. 
The  actual  crustal  thickness  in  the  region  of  th  south  end  is  about  25  km  [  Braile  et  al., 
1989  ],  10  km  thinner  than  that  was  used.  Thus,  a 1  or  /?  ’  is  about  one  second  greater 
than  a  and  /?,  the  travel  time  is  smaller  for  any  ray  ending  within  this  region. 

Error  analysis  of  azimuthal  anisotropy 

The  azimuthal  velocity  residuals  plotted  in  Figure  9  are  calculated  the  Pn  velo¬ 
city  image  given  in  Figure  7,  without  correcting  the  mantle  velocity  gradient.  The 
averaged  values  are  almost  identical  to  that  derived  from  a  model  with  a  positive  velo¬ 
city  gradient  mantle.  This  implies  that  the  assumption  of  a  homogeneous  mantle  does 
not  effect  the  velocity  residual  dependence  on  azimuth.  However,  the  two  station 
method  of  Beghoul  and  Barazangi  [1990]  tends  to  enlarge  the  contribution  of  the  man¬ 
tle  velocity  gradient  to  the  velocity  residuals.  For  example,  the  Pn  velocities  from  the 
two  station  method  are  7.82  km/fe  for  distance  of  500  km  (  two  stations  at  1000  and 
1500  km  )  and  7.73  kmyfe  for  distance  of  300  km  (  two  stations  at  200  and  500  km) 
instead  of  the  real  velocity  of  7.72  (200  km)  and  7.73(500  km),  using  v0  7.72  km^ 
and  c  2.56  XlO^s-*  obtained  previously.  The  correction  of  the  average  mantle  velo¬ 
city  gradient  should  not  affect  the  azimuthal  anisotropy,  since  the  correction  is  the 
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same  in  every  direction. 

Assuming  a  homogeneous  model  for  both  crust  and  mantle,  we  obtained  the 
velocity  image  given  in  Figure  8,  the  azimuthal  velocity  residuals  of  which  are  given 
in  Figure  10.  The  similarity  of  Figure  5  to  Figure  10  infers  that  the  initial  crustal 
models  do  not  effect  the  azimuthal  dependence  of  the  velocity  residuals.  Note  that  the 
contributions  of  the  lateral  variation  of  the  Pn  velocities,  the  mislocation  of  the  events, 
and  lateral  variation  of  the  crustal  structure  to  the  azimuthal  Pn  velocity  residuals  were 
taken  out  in  the  course  of  inversion  and  allocation  procedure. 

We  use  the  refined  Pn  travel  time  data  set,  which  is  not  yet  put  through  the  last 
inversion  to  get  the  tomographic  image  in  Figure  8.  This  means  that  the  lateral  varia¬ 
tions  in  the  Pn  velocities  is  put  back  to  the  travel  times,  and  the  a  and  /?  values  are 
refined.  The  average  velocity  of  7.77  ±0.16  km/fe  was  given  above  for  this  data  set  ( 
Figure  8  ).  We  use  this  values  as  the  velocity  for  the  upper  mantle  beneath  the  whole 
studied  area  (  homogeneous  ).  The  azimuthal  velocity  residuals  are  given  in  Figure  11. 
The  fluctuation  of  the  averaged  velocity  residuals  is  larger  than  the  residuals  shown  in 
Figure  5. 

Figure  12  gives  the  azimuthal  velocity  residuals  for  the  raw  data  with  correct¬ 
ing  a  homogeneous  crust  and  upper  mantle.  The  average  velocity  is  7.77  ±0.20 
kmyhec.  The  standard  error  is  larger  than  that  in  Figure  11.  The  pattern  of  Figure  12  is 
roughly  the  same  as  that  of  Figure  11.  This  figure  suggests  that  as  long  as  the  ray 
paths  are  dense  there  is  no  evidence  for  azimuthal  anisotropy  in  the  Basin  and  Range 
province,  regardless  of  the  assumptions  made  in  the  crustal  and  mantle  velocity  struc¬ 
tures. 


Discussion  and  conclusions 


Pn  velocity  image 


-  14- 


The  inversion  method  and  procedure  used  in  this  study  to  get  the  Pn  velocity 
image  differ  from  those  in  previous  Pn  tomography  studies,  for  example,  Hearn  et  al. 
[1901].  We  used  the  following  two  steps  prior  to  the  back  projection  inversion: 

(1)  We  estimated  the  vertical  travel  times  for  the  sources  and  receivers,  with  a  priori 

knowledge  on  the  lateral  variation  of  in  crustal  structure.  These  delay  times  were 
then  subtracted  from  the  Pn  travel  times  to  obtain  the  apparent  Pn  velocity. 

(2)  We  estimated  the  average  mantle  velocity  gradient,  which  was  used  to  correct  the 

apparent  P„  velocities  to  get  P  velocities  along  the  uppermost  mantle. 

The  numerical  experiment  and  theoretical  development  that  make  these  two 
corrections  possible  are  given  in  this  study. 

The  average  Pn  velocity  from  the  inversion  is  7.72  ±0.16  kmyfeec.  The  average 
mantle  velocity  gradient  obtained  in  this  study  is  2.0  XIO-V^,  including  the  contribu¬ 
tion  of  the  Earth’s  sphericity,  1.2  XlO^s-1.  Two  distinct  low  Pn  velocity  regions 
beneath  the  Sierra  Nevada  mountains  and  southern  Nevada  seem  to  be  connected  and 
form  one  large  low  velocity  region  in  the  southern  Basin  and  Range  province.  The  Pn 
velocity  is  less  than  7.8  kmyfe  in  the  south;  and  it  is  greater  than  the  7.8  kniyfe  in  the 
north.  These  slow  and  fast  regions  form  two  big  eastrwest  strips  (  Figure  4  ).  The  Pn 
velocity  pattern  of  low-south  and  f as tr north,  in  an  almost  east-west  in  direction,  could 
be  due  to  upwelling  beneath  the  Sierra  Nevada  and  crustal  extension  in  the  southern 
Nevada,  and  as  discussed  above.  It  may,  also,  represent  a  small  convection  current 
occurred  here  in  the  region,  upwelling  in  the  south  and  downwelling  in  the  north 
beneath  the  Basin  and  Range  province,  if  the  velocity  differences  are  caused  by  the 
temperature. 

The  Pn  velocity  image  is  not  affected  much  by  the  mantle  velocity  gradient 
But,  the  three  dimensional  crustal  structure  is  crucial  to  the  Pn  velocity  image.  The 
lower  or  faster  velocity  pattern  is  largely  the  result  of  the  too  short  or  too  long  esti¬ 
mates  of  the  station  and  receiver  delay  times.  The  resolution  of  the  velocity  image  is 


-  15  - 


generally  within  three  degrees,  within  two  degrees  for  the  grids  with  denser  ray  cover¬ 
age,  and  within  five  degrees  for  poorer  ray  coverage  grids. 

Azimuthal  anisotropy  of  P&  velocities 

The  method  used  in  this  study  to  infer  the  Pn  velocity  anisotropy  also  differs 
from  previous  studies,  e.  g.,  Vetter  and  Minster  [  1081  ],  in  the  following  ways: 

(1)  Most  important  of  all,  we  took  away  the  lateral  variation  in  PD  velocities  from  the 

possible  azimuthal  dependence  of  the  velocity  residuals,  by  using  tomographic 
inversion. 

(2)  Relocation  of  the  sources  and  refinement  of  the  local  crustal  structures  beneath  the 

sources  and  receivers. 

Pn  azimuthal  anisotropy  for  this  region,  if  any,  is  not  likely  to  be  greater  than 
one  percent  from  velocity  residuals  calculated  from  various  models.  However,  previ¬ 
ous  Pn  studies  suggested  more  than  three  percent  anisotropy  [  Vetter  and  Minster, 
1981;  Beghoul  and  Barazangi,  1990  ].  The  obvious  reason  for  this  discrepancy  is  due 
to  the  heterogeneity  and  the  much  larger  volume  of  the  present  dataset  An  azimuthal 
dependence  might  show  up  if  the  stations  were  located  within  some  limited  region  as 
the  case  considered  by  Vetter  and  Minster  [  1981  ].  For  example,  all  stations  are 
located  in  the  center  of  the  region  with  earthquakes  around  them  (  Figure  8).  The 
seismic  anisotropy  in  this  method  is  merely  the  lateral  variation.  The  two  station 
method,  which  was  used  by  Beghoul  and  Barazangi  [  1990  ]  to  estimate  P&  velocities 
for  this  region,  does  not  account  for  the  mantle  velocity  gradient.  Mantle  velocity  gra¬ 
dient  and  lateral  variations  of  Pn  velocities  are  responsible  for  the  anisotropy  in  this 
case. 

For  the  azimuthal  velocity  residuals,  the  effects  of  the  mantle  velocity  gradient 
and  the  initial  crustal  models  are  small.  The  lateral  variation  of  the  Pn  velocities  and 
refinements  of  the  vertical  travel  times  of  the  sources  and  receivers  are  substantial; 


-  16- 


they  affect  both  the  amplitude  and  general  pattern. 

Savage  et  al.,  [  1690  ]  found  favorable  evidence  of  azimuthal  anisotropy  from 
shear  wave  splitting  studies.  This  may  suggest  that  there  is  no  anisotropy  a  few  tens  of 
kilometers  beneath  Moho,  but  at  greater  depths  anisotropy  appears,  since  the  ray  paths 
for  Pn  are  horizontal  while  that  of  SKS  are  almost  vertical.  It  may  also  suggest  that  the 
mineralogy  at  the  shallow  depth  beneath  Moho  has  a  rich  content  of  some  mineral 
which  has  a  strong  shear  velocity  anisotropy  and  a  weak  compressional  velocity  aniso¬ 
tropy,  such  as  spinel  and  rutile. 

In  conclusion,  we  applied  the  tomographic  inversion  method  to  1226  Pn  travel 
times  from  ISC  catalogs  to  study  the  Pn  velocity  distribution  and  anisotropy  beneath 
the  Basin  and  Range  province.  We  obtained  a  velocity  image  with  slow  velocities  in 
the  south  and  fast  velocities  in  the  north.  The  image  may  suggest  a  upwelling  beneath 
the  Sierra  Nevada  and  a  crustal  extension  beneath  southern  Nevada.  It  may  also  sug¬ 
gest  a  small  convection  current  occurred  beneath  this  region,  with  upwelling  in  the 
south  and  downwelling  in  the  north.  Before  the  inversion,  we  refined  the  apparent  Pn 
velocities  by  correcting  vertical  travel  times  for  sources  and  receivers  using  three 
dimensional  crustal  structures  and  correcting  for  the  mantle  velocity  gradient,  which 
are  progressively  refined  in  the  couree  of  the  inversion  and  allocation  procedure.  The 
aforementioned  refining  process  of  the  dataset  is  crucial  in  obtaining  an  accurate  Pn 
velocity  image,  but  does  not  affect  the  azimuthal  dependence  of  the  velocity  residuals, 
although  the  amplitude  of  these  is  increased.  Azimuthal  anisotropy  in  the  uppermost 
mantle  beneath  the  Basin  and  Range  province  is  veiy  small,  less  than  one  percent,  if 
present  at  all,  in  spite  of  the  strong  orientation  of  the  local  stress  fields. 
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Appendix 


We  assume  that  the  mantle  velocity  is  given  by  Equation  (3),  i.  e., 

v(z)  -  t>0(l  +  cz)  .  (Al) 

For  head  waves  (  such  as  Pn)  turning  at  a  depth,  h,  the  ray  parameter  is  1  /v(h), 
and  the  travel  time  T(h )  is  given  by: 


T(h) 


v(h)  o  v(z)2  v(h)2 

-  V0(i+ct )  [X~7«1+a );-l)l/2+'-(1^'1 }  ■■>(!+<*  +((l+cA )2-l )»/«)]  ,(A2) 

where  X  is  the  horizontal  distance  traveled  in  the  mantle.  The  head  wave  travel 
time  over  X  should  satisfy  dT(h)/dh  «  0,  we  have: 


X  =  ±({l+ch)2-l)'/2  , 
c 

and  the  head  wave  travel  time  over  X  is: 


(A3) 


VnC 


ln(l+cA  +  ((1+cA  )2— l)1/2) . 


(A4) 


From  definition,  The  observed  head  wave  velocity  is  vobe  «=X/!T  .  From 
(A3)  and  (A4)  we  have: 

_  y/2ch  (1+ch /4)  _  V2ch  (1+ch /4)  _  ,,  ,  cL  /4r, 

0  Tn (l+ek  +V2dT  (/+cA/4))  ~V°  /12)  ~t'o(1+T)  ’  (A5) 

if  ch  «1.  For  ch  «1,  X  can  be  approximated  by 

X  ~  V8 h /c  ,  or  h  ~  cX2/ 8  .  (A6) 

Substituting  (A6)  into  (A5),  we  have  Eq.(4)  in  the  text. 

vobe  w  vo(!  +  C2X2/ 24)  . 

If  X>Xm~((l+c#r2)—  l)1//2/c,  at  which  the  ray  bottoms  the  depth  of  Hm, 
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below  which  Eq.(Al)  does  not  hold,  Eq.  (A4)  should  be  written  as: 

«o(l+ctfm)  «0C 

_  X-Xm  1  x„ 

+  Vo(1  +  icffm) 

_  X  +  ^XmcHm 
*  ^l  +  «ff.)  ' 

From  the  definition  of  voia  «■  X/T,v/e  have 


vo6«  * 


2Xm 

l  +  ~sTcH” 

v0{1  +  cHm ) 


ti  ft  \ 

«  v0  +  v0ctfm(l — , 


(A7) 


Eq.(5)  in  the  text. 
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Figure  1.  Ray  paths.  Triangles  are  stations  and  stars  are  earthquakes.  Heavy  dash  lines 
depict  the  boundaries  of  the  Basin  and  Range  province.  State  boundaries  are  also 
included. 


Figure  2.  Two  dimensional  model  used  to  generate  finite  difference  synthetic  seismo¬ 
grams.  Ibe  velocities  for  blocks  are  given  in  the  boxes.  The  velocity  is  5.2 
kmyfeec  for  the  layer  at  the  top,  for  the  layer  beneath  it,  the  velocity  is  6.0 
kmyfeec.  The  star  is  the  source. 


Figure  3.  Comparison  of  the  theoretical  calculations  using  Eqs.  (1)  and  (2)  with  the 
numerical  calculation  using  finite  difference  synthetic  seismogram  method.  The 
solid  lines  are  truncated  synthetic  seismograms.  The  solid  dots  are  the  travel  time 
picks  from  the  truncated  seismograms,  and  the  arrows  are  the  theoretical  calcula¬ 
tions.  The  model  used  is  given  in  Figure  2. 


Figure  4.  Pn  velocity  image  for  the  Basin  and  Range  province.  Heavy  dash  lines  depict 
the  boundaries  of  the  Basin  and  Range  province.  The  diagonal  ruled  area  are  the 


area  without  enough  resolution. 


Figure  5.  Velocity  residuals  (  6v/v  in  percent )  calculated  by  using  Equation  (11), 

and  average  velocity  residuals  (  Equation  (12)  )  in  the  20  degree  interval  (  solid 
dot )  with  error  bars  (  square  root  of  Equation  (13) ) ,  versus  the  azimuths  of  the 


Figure  6.  Point  spread  functions  (  p-s.f. )  for  six  grids,  five  close  to  the  boundaries,  one 
in  the  center  of  the  region.  Without  losing  any  information,  I  put  three  of  them 
on  one  plot,  since  they  do  not  overlap.  p.s.f.  of  a)  middle  grid  2  degrees,  left 
grid  5  degrees,  and  right  grid  3  degrees;  and  b)  left  grid  spreads  2  degrees,  lower 
grid  2  degrees,  and  right  grid  2.5  degrees.  Heavy  dash  lines  depict  the  boundaries 
of  the  Basin  and  Range  province.  The  diagonal  ruled  area  are  the  area  without 
enough  resolution,  or  with  the  amplitude  less  than  0.1. 


Figure  7.  Pn  velocity  image  for  the  Basin  and  Range  province  without  correcting  for 
mantle  velocity  gradient*  i.  e.,  assuming  7'  “0.0  seconds.  The  diagonal  ruled 
area  are  the  area  without  enough  resolution. 


Figure  8.  Pn  velocity  image  for  the  Basin  and  Range  province  assuming  a  homogene¬ 


ous  crust  and  mantle. 
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Figure  0.  Velocity  residuals,  averaged  velocity  residuals  after  correcting  the  velocity 
distribution  of  Figure  7,  versus  the  azimuths  of  the  paths. 
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Figure  10.  Velocity  residuals,  averaged  velocity  residuals  after  correcting  the  velocity 
distribution  of  Figure  8,  versus  the  azimuths  of  the  paths. 
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Figure  11.  Velocity  residuals,  averaged  velocity  residuals  before  correcting  the  velocity 
distribution  of  Figure  8,  versus  the  azimuths  of  the  paths. 
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Lateral  Variation  in  Crustal  Structure  of  the  Northern  Tibetan 
Plateau  Inferred  from  Teleseismic  Receiver  Functions 

by  Lupei  Zhu,  Thomas  J.  Owens,  and  George  E.  Randall 


Abstract  A  one-year  period  of  passive-source  recordings  on  the  Tibetan  Plateau 
includes  more  than  300  teleseismic  events  ranging  from  35°  to  95°  in  epicentral  dis¬ 
tance  and  with  good  back-azimuth  coverage.  This  provides  a  chance  to  study  the 
lateral  variation  of  the  crustal  structure  along  the  northern  boundary  of  the  plateau  by 
investigating  the  variation  of  source-equalized  P  waveforms  (receiver  functions)  with 
back-azimuth  at  stations  BUDO,  ERDO  and  TUNL  located  in  this  region.  It  is  shown 
clearly  that  the  waveforms  of  first  5  sec  receiver  functions  vary  with  back-azimuth  sys¬ 
tematically:  the  radial  components  are  symmetric  across  N-S  axis  while  the  tangential 
components  are  mirror-symmetric  across  this  axis.  All  these  can  be  modeled  well  by 
E-W  striking  dipping  interfaces  in  the  upper-middle  crust.  The  strike  direction  is  very 
consistent  with  the  E-W  trend  of  surface  geology.  Modeling  a  P-to-S  converted  phase 
in  the  radial  components  at  each  station  shows  that  there  is  a  mid-crust  low  velocity 
layer  with  its  upper  interface  dipping  20-30°  to  the  south.  In  addition,  a  shallow 
northwards-dipping  interface  under  ERDO  is  responsible  for  the  "double  peaked"  direct 
P  arrivals  in  the  radial  components  from  certain  back-azimuths  and  large  tangential 
components.  The  low  velocity  layer,  together  with  other  geology  and  seismology 
observations,  suggests  that  there  is  a  hot,  possible  partial  melting  zone  in  the  middle 
crust  under  the  northern  plateau.  Alternately  dipping  velocity  interfaces  might  be  asso¬ 
ciated  with  some  buried  thrust  faults  in  the  brittle  upper  crust  which  accommodated 
crust  shortening  during  some  early  stage  of  the  plateau  formation. 

Introduction 

The  Tibetan  Plateau,  bounded  by  the  Kunlun  mountains  in  the  north  and 
Himalayas  in  the  south,  with  its  outstanding  average  elevation  of  5  km  and  an  excep¬ 
tional  area  extent  of  nearly  7xl05km2,  has  long  been  an  interesting  and  challenging 
subject  in  geoscience.  Although  it  is  now  generally  agreed  that  the  uplift  of  Tibet  pla¬ 
teau  is  the  result  of  the  collision  of  Indian  plate  with  Asia  during  Middle  Eocene  45 
million  years  ago,  after  the  Tethys  ocean  lithosphere  subducted  completely  beneath  the 
Asia,  the  details  of  this  collision  mechanism  and  the  plateau  uplift  still  remain  contr¬ 
oversial.  Several  collision  models  have  been  proposed.  Among  the  most  influential 
ones  are  the  underthrust  model  and  its  variants  (Argand,  1924;  Zhao  and  Morgan, 
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1985;  Beghoul  et  al.,  1993),  the  crust  shortening  and  thickening  model  (Dewey  and 
Bird,  1970;  Dewey  and  Burke,  1973),  and  the  lateral  crustal  extrusion  model  (Molnar 
and  Tapponmer,  1977;  Tapponnier  et  al.,  1982).  Due  to  lack  of  necessary  constraints 
on  the  lithospheric  structure  under  the  plateau,  it  is  difficult  to  discriminate  among 
different  models. 

So  far,  most  geological  and  geophysical  investigations  on  the  plateau  were  con¬ 
centrated  on  the  southern  Tibet  near  the  Himalayas.  Comparatively,  the  structure  near 
the  northern  boundary  of  the  plateau  was  poorly  resolved.  Since  the  northern  part  of 
the  plateau  should  play  an  equally  important  role  as  its  southern  counterpart  during  the 
forming  of  the  plateau,  better  understanding  of  the  crust  and  upper  mantle  structure 
under  the  northern  plateau  will  help  to  constrain  dynamic  models  of  plateau  formation. 

Recently,  as  a  joint  research  project  between  Institute  of  Geophysics,  State 
Seismological  Bureau,  China,  and  University  of  South  Carolina,  State  University  of 
New  York  at  Binghamton,  1 1  broadband  3-component  seismic  recorders  were 
deployed  on  the  plateau  for  one  year  of  recording  (Fig.  1).  Three  stations  (TUNL, 
BUDO  and  ERDO)  were  placed  across  the  northern  boundary  of  the  plateau,  providing 
an  ideal  opportunity  to  study  the  crust  and  upper  mantle  structure  in  this  area.  In  this 
paper,  we  will  use  teleseismic  P  waveform  data  to  study  the  lateral  variation  of  the 
crust  structure.  We  will  show  that  there  is  strong  lateral  variation  in  the  crust  of 
northern  Tibet  and  that  most  of  this  variation  is  caused  by  E-W  striking  dipping  velo¬ 
city  interfaces  in  the  upper-middle  crust. 

Method 

Teleseismic  P  waveform  data  contains  information  about  both  the  source  side  and 
the  receiver  side  structures.  Using  the  method  proposed  by  Langston  (1979),  the  source 
and  path  effects  can  be  removed  from  the  P  waveforms  by  deconvolving  the  vertical 
component  from  the  horizontal  components  under  the  assumption  that  the  vertical 
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component  of  the  teleseismic  recording  consists  primarily  of  the  unwanted  source  and 
path  effects.  The  source-equalized  teleseismic  waveform  which  is  called  the  receiver 
function  is  most  sensitive  to  the  P-to-S  conversions  from  structures  beneath  the  record¬ 
ing  site  (Owens  et  al.,  1984). 

The  radial  component  of  broadband  receiver  functions  is  often  modeled  by  hor¬ 
izontally  stratified  velocity  structure  using  a  time-domain  inversion  technique  (Owens 
et  al.,  1984;  Ammon  et  al.,  1990).  This  method  has  been  applied  to  numerous  areas 
and  proven  to  be  efficient  in  obtaining  the  velocity-depth  variation  of  the  crust  and  the 
Moho  depth.  However,  if  the  tangential  component  of  receiver  function  is  comparable 
in  amplitude  to  the  radial  component,  which  generally  indicates  that  there  is  strong 
lateral  heterogeneity  within  the  crust  beneath  the  site,  modeling  the  details  of  the  radial 
component  with  horizontally  stratified  velocity  will  lead  to  erroneous  results. 

One  of  the  special  cases  of  lateral  heterogeneity  is  the  dipping  planar  velocity 
interface  in  the  crust.  In  this  case,  both  the  radial  and  tangential  component  vary  with 
back-azimuth  in  a  quite  predictable  pattern.  Figure  2  shows  the  variation  of  synthetic 
receiver  functions  for  a  crust  model  which  contains  a  mid-crustal  low  velocity  layer 
with  its  upper  interface  dipping  20°  to  the  south.  The  interface  is  at  depth  of  20  km 
where  the  S  velocity  jumps  from  3.0  km/sec  below  it  to  3.5  km/sec  above.  The  Moho 
is  at  a  depth  of  65  km.  Teleseismic  ray  parameter  is  0.068  sec/km  which  corresponds 
epicentral  distance  of  50°.  Synthetics  are  calculated  by  using  a  3D  ray-tracing  method 
of  Langston  (1977).  It  can  be  seen  that  the  waveform  variation  with  back-azimuth  is 
very  systematic:  the  radial  components  are  symmetric  across  a  line  parallel  to  the  dip 
direction  of  the  interface,  while  the  tangential  components  are  mirror-symmetric  about 
this  direction  (indicated  by  arrow  in  the  figure).  Tangential  component  has  the  largest 
energy  in  the  strike  direction  and  will  be  zero  when  ray  comes  from  direction  parallel 
to  the  dip  direction.  The  P-to-S  converted  wave  generated  at  the  dipping  interface  is 
labeled  as  phase  A.  This  negative-amplitude  phase  emerges  in  the  southern  back- 
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azimuths  with  maximum  amplitude  from  the  dip  direction  and  vanishes  gradually  away 
from  the  dip  direction. 

So,  if  the  dipping  interface  is  the  major  cause  of  lateral  heterogeneity  under  the 
station,  it  is  very  easy  to  determine  the  dip  and  strike  direction  by  looking  at  the 
back-azimuthal  profile  of  the  receiver  functions.  The  dip  angle,  velocity  contrast,  and 
the  depth  of  the  interface  can  be  obtained  by  modeling  the  amplitude  and  arrival  time 
of  the  converted  phase  from  the  interface,  although  a  trade-off  exists  between  the  dip 
angle  and  the  velocity  jump  across  the  interface  (Owens  et  al.,  1988).  Figure  3a  shows 
the  variation  of  absolute  amplitude  of  the  P-to-S  converted  phases  from  above  dipping 
interface  with  different  dip  angles.  Figure  3b  is  the  variation  of  maximum  amplitude 
with  velocity  contrast  at  the  interface.  It  can  be  seen  that  the  dip  angle  controls  the 
shape  of  the  amplitude  variation  with  back  azimuth,  while  the  velocity  difference 
across  the  interface  contributes  mainly  to  the  maximum  amplitude  (Fig.  3).  The  dip 
angle  also  affects  the  maximum  amplitude.  But  to  produce  same  maximum  amplitude 
with  velocity  contrast  differing  0.1  km/sec,  the  dip  angles  will  differ  more  than  10° 
(Fig.  3b),  which  will  have  quite  different  shapes  of  amplitude  variation  with  back- 
azimuth  (Fig.  3a).  For  this  reason,  if  data  set  has  enough  back-azimuth  sampling,  rea¬ 
sonable  constraints  can  be  made  on  the  dip  angle  and  velocity  contrast  of  the  interface. 

Data 

Ten  of  the  11  stations  in  the  experiment  were  equipped  with  Streckeisen  STS-2 
sensors.  Station  TUNL  was  equipped  with  Guralp  CMG-3ESP  sensor.  The  STS-2  has  a 
broad  velocity  response  with  comer  frequencies  of  0.0083  Hz  and  50  Hz  and  the 
CMG-3ESP  has  comer  frequencies  of  0.033  Hz  and  30  Hz.  40  samples/sec  sampling 
rate  was  used  for  recording  waveforms  of  regional  and  teleseismic  events.  Details 
about  the  instrumentation  of  the  experiment  can  be  found  in  Owens  et  al.,  1993. 
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More  than  300  teleseismic  events  ranging  between  35°  to  95°  in  epicentral  dis¬ 
tance  were  recorded  during  the  Tibetan  experiment.  A  preliminary  study  of  the 
receiver  functions  of  teleseismic  events  from  certain  back-azimuths  on  3  stations 
(TUNL,  WNDO  and  XIGA)  located  in  the  northern,  central  and  southern  plateau  has 
shown  the  great  thickness  of  crust  (60  to  80  km)  and  the  difference  of  velocity  struc¬ 
tures  between  north-central  and  southern  Tibet  (Zhu  et  al.,  1992).  Large  tangential 
components  were  also  observed  in  TUNL  and  XIGA. 

Three  stations,  TUNL,  BUDO,  and  ERDO,  were  located  near  the  north  boundary 
of  the  plateau  which  is  defined  by  an  E-W  striking  fault  system,  the  Kunlun  Mountain 
Fault  System  (Fig.  1).  The  elevation  increases  from  3  km  in  the  Qaidam  Basin  to 
more  than  4.5  km  inside  the  plateau.  The  surface  geology  in  this  region  includes 
highly  deformed  Permian  subduction-related  rocks  and  Upper  Paleozoic  pluton  rocks 
(Dewey  et  al.,  1988).  The  Kunlun  fault  system  consists  of  several  large  left-lateral 
strike-slip  faults.  Some  thrust  faults  were  discovered  during  several  recent  geotraverses, 
for  example,  the  South  Qaidam  Border  Thrust  Fault  and  Fenghuoshan  thrust  faults 
(Molnar  et  al.,  1987;  Kidd  et  al.,  1988). 

Figure  4  is  the  distribution  of  teleseismic  events  used  in  this  study.  The  good 
back-azimuthal  coverage  of  teleseismic  events  enables  us  to  study  the  lateral  variation 
of  the  crustal  structure  in  the  area  from  the  variation  of  the  receiver  function 
waveforms  with  back-azimuth  at  each  of  these  stations.  The  3 -component  teleseismic 
recordings  were  cut  with  a  time  window  of  180  sec  length  beginning  60  sec  before  P 
arrival.  Then  horizontal  components  were  rotated  to  radial  and  tangential  directions 
and  deconvolved  with  vertical  component  in  frequency  domain.  Gaussian  parameter  of 
2.5  was  used  to  exclude  frequencies  above  1  Hz.  To  make  the  deconvolution  algorithm 
stable,  we  used  a  range  of  water-levels  from  0.1  to  0.0001  and  chosed  the  best  one 
from  them.  Details  on  the  calculation  of  receiver  function  are  described  in  some  pub¬ 
lished  papers  (Owens  et  al.  1984;  Ammon  et  al,  1990). 
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For  each  station,  the  teleseismic  events  are  divided  into  groups  of  back-azimuth 
range  less  than  10°,  according  to  their  natural  distribution  with  respect  to  that  station 
and  the  waveform  similarity  of  the  receiver  functions.  We  found  that  the  difference  in 
epicentral  distance  does  not  affect  the  waveforms  much,  so  some  of  the  groups  span  a 
distance  range  of  about  20°.  The  receiver  functions  in  each  group  are  then  stacked. 
The  back-azimuthal/distant  range  and  the  number  of  events  for  each  stacking  suite  are 
listed  in  Table  1. 

Figure  5  is  the  back-azimuthal  profiles  of  the  stacked  receiver  functions.  It  is 
clear  that  there  is  strong  lateral  variation  in  the  upper-middle  crustal  structure  of  this 
region  as  evidenced  by  the  large  tangential  components  in  the  first  5  sec  on  all  these  3 
stations.  Also  the  waveforms  vary  with  back-azimuth  in  the  pattern  predicated  by  dip¬ 
ping  interfaces  we  described  above.  Both  the  symmetry  of  radial  components  and  the 
mirror-symmetry  of  tangential  components  across  N-S  direction  are  obvious  in  the  first 
5  sec  of  the  receiver  functions.  The  tangentials  are  nearly  zero  in  the  N-S  direction  (as 
indicated  by  arrows  in  Figure  5).  Largest  tangential  components  are  on  events  arriving 
from  the  east  and  west  directions  (back  azimuths  of  about  90-100°  and  270-280°)  with 
opposite  polarities. 


Modeling  Results 

One  prominent  feature  observed  in  the  radial  receiver  functions  is  a  mid-crust  P- 
to-S  converted  phase  in  back-azimuth  range  of  about  100°  to  300°  (indicated  by 
dashed-lines  in  Figure  5).  This  negative  polarity  arrival  is  well  illustrated  at  BUDO 
(time  delay  3.2  sec  with  respect  to  direct  P  arrival)  and  also  can  be  identified  at  TUNL 
(2.4  sec)  and  ERDO  (3.8  sec).  Since  there  is  no  large  converted  phase  before  it,  this 
phase  could  not  be  multiple  converted  phase.  So  its  negative  polarity  indicates  that  it  is 
generated  at  top  of  a  lower  velocity  layer  (as  shown  in  synthetic  receiver  functions  in 
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Fig.  2).  The  variation  of  the  amplitude  of  this  phase  with  back-azimuth  is  caused  by 
the  southward  dipping  of  the  interface. 

We  measure  the  absolute  amplitudes  of  this  phase  at  the  radial  components  from 
different  back-azimuths  and  fit  them  with  synthetical  calculations  by  trial-and-error  for¬ 
ward  modeling.  A  range  of  S-velocity  contrasts  from  3.4/3. 1  to  3.6/2.8  across  the 
interface  and  dip  angles  from  10°  to  30°  have  been  tested.  We  use  maximum  ampli¬ 
tude  to  determine  the  velocity  contrast  and  dip  direction  and  use  the  shape  of  the 
amplitude  variation  with  back-azimuth  to  determine  the  dip  angle.  Figure  6  gives  the 
modeling  results.  There  is  very  good  back-azimuth  sampling  at  BUDO  which  leads  to 
substantial  constraint  on  the  velocity  contrast  (3 .6/2.9)  and  dip  angle  (25°)  (Fig.  6b). 
Station  TUNL  and  ERDO  do  not  have  enough  data  points  to  resolve  the  dip  angles. 
The  velocity  contrast  is  3.5/3.0  across  the  interface  and  the  dip  angle  can  be  25±5°  (at 
TUNL)  or  20±5°  (at  ERDO).  The  depth  of  the  mid-crustal  dipping  interface  at  each 
station  is  determined  by  the  arrival  time  of  the  converted  phase,  assuming  the  average 
upper  crust  S-velocity  of  3.5  km/sec.  The  results  are  15  km  for  TUNL;  24  km  for 
BUDO  and  30  km  for  ERDO. 

Another  feature  of  radial  receiver  functions  is  the  apparent  time  shift  of  direct  P 
arrivals  from  back-azimuths  30°  and  315°  at  ERDO  (indicated  by  blacken  dots  in  Fig. 
5c).  It  is  actually  produced  by  a  shallow  dipping  interface  which  causes  the  direct  P 
arrival  very  small  and  the  Ps-wave  generated  on  the  interface  appears  to  be  the  first 
arrival.  Such  shallow  dipping  interface  effect  on  the  receiver  function  has  been 
modeled  and  discussed  by  Owens  and  Crosson  (1988).  We  found,  by  adjusting  the 
depth,  dip  angle  and  velocity  contrast  of  the  interface,  that  a  low  velocity  top  layer 
(S-velocity  2.5  km/sec)  of  thickness  4  km  with  bottom  interface  dipping  15°  to  the 
north  can  model  both  the  radial  and  tangential  waveform  of  first  2-3  sec  very  well 
(Fig.  7) 
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Conclusions  and  Discussion 

In  summary,  the  waveform  variation  of  receiver  functions  with  back-azimuth  at 
three  stations  located  in  the  northern  Tibetan  Plateau  shows  that  there  are  strong  lateral 
heterogeneities  within  the  upper-middle  crust  in  this  region.  The  waveform  variation 
pattern  also  indicates  that  it  is  mainly  produced  by  E-W  striking  dipping  interfaces. 
The  strike  direction  is  very  consistent  with  the  E-W  trend  of  surface  geology.  Model¬ 
ing  a  negative-polarity  P-to-S  converted  phase  in  the  radial  component  at  each  station 
shows  that  there  is  a  mid-crust  low  velocity  layer  with  upper  interface  dipping  20-30° 
to  the  south.  In  addition,  a  shallow  northwards-dipping  interface  under  ERDO  is 
responsible  for  the  "double  peaked"  direct  P  arrivals  in  the  radial  components  from 
certain  back-azimuths  and  large  tangential  components. 

Our  result  of  mid-crust  low  velocity  layer  (LVL)  is  consistent  with  previous 
results  from  phase  velocity  dispersions  of  surface  wave  (Romanowicz,  1982;  Brandon 
and  Romanowicz,  1986).  The  modeling  shows  that  the  shear  velocity  in  the  middle 
crust  could  by  as  low  as  2.9  km/s.  Due  to  difficulty  of  identifying  converted  phase 
from  bottom  of  the  LVL,  we  are  unable  to  determine  the  thickness  of  the  LVL.  But 
the  LVL  itself  has  important  indication  on  the  crustal  structure  of  the  plateau.  One  of 
possibilities  to  get  such  low  shear  velocity  at  depth  larger  than  20  km  is  partial  melt¬ 
ing.  Partial  melting  is  effective  in  decreasing  velocities  under  high  presure  (Anderson, 
1989).  The  very  lower  average  S-velocity  (less  than  3.5  km/s)  from  surface  wave 
dispersion  data  led  Chen  and  Molnar  to  suggest  that  lower  crust  of  the  northern  plateau 
is  in  partial  melting  (Chen  and  Molnar,  1981).  Recent  earthquake  relocations  on  Tibet 
plateau  show  that  most  earthquakes  occur  only  in  the  depth  range  between  5  and  15 
km  (Chen  and  Molnar,  1983;  Zhao  and  Helmberger,  1991).  The  partial  melting  is  also 
supported  by  the  observation  of  Quaternary  volcanic  activities  in  the  northern  Tibet 
(Deng,  1978;  Pearce  and  Mei,  1988).  The  high  temperature  in  the  crust  could  be  due 
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to  the  concentration  of  crust  material  during  the  crust-shortening  or  high  heat  flow 
from  the  upper  mantle. 

The  E-W  striking  dipping  interfaces  of  our  modeling  results  might  be  associated 
with  some  buried  thrust  faults  in  the  brittle  upper  crust  in  this  northern  margin  of  the 
plateau.  For  example,  the  shallow  northwards-dipping  interface  under  ERDO  could  be 
related  to  the  Fenghuoshan  thrust  faults  which  run  EW  in  the  south  of  the  station.  It  is 
estimated  that  large  amount  of  crust  shortening  has  occurred  during  the  formation  of 
the  plateau.  Palaeomagnetic  data  provided  a  rough  relative  movement  between  the 
India  plate  and  the  ’stable’  Eurasia.  It  shows  that  only  small  fraction  of  India  crust 
(400  km)  was  absorbed  since  the  collision  and  most  of  this  fraction  can  be  found  accu¬ 
mulated  on  the  Himalayan  front  instead  of  underthrusting  into  the  Tibet  crust  (Dewey 
et  al.,  1988).  On  the  other  hand,  about  2000  km  crust  shortening  has  occurred 
between  the  Indo-Zanbo  collision  suture  and  the  Siberia  (Patriat  and  Achache,  1984; 
Lin,  1988).  Accommodating  such  large  amount  of  crust  shortening  needs  thrust  faults 
within  the  plateau  or  along  the  north  and  south  boundaries.  Large  scale  thrust  faults 
exist  along  the  Himalayan  front  and  north-west  boundary  between  Tibet  and  Tarim 
Basin.  Recent  geological  expeditions  also  found  evidence  of  thrust  faults  in  northern 
Tibet  (Molnar  et  al.,  1987;  Tapponnier  et  al.,  1990). 

Shallow  and  mid-crust  dipping  interfaces  have  noticeable  impact  on  the  modeling 
of  deep  structure  such  as  the  Moho.  As  demonstrated  in  our  synthetic  receiver  func¬ 
tions  in  Figure  2,  the  multiple  converted  phases  of  the  dipping  interface  will  interfere 
with  the  Moho  Ps  and  make  the  amplitude  of  Moho  Ps  varying  with  back  azimuth. 
This  effect  is  also  seen  in  the  observed  receiver  functions.  For  crust  structure  of  Tibet 
with  thickness  of  about  65  km,  the  predicted  arrival  time  of  Moho  Ps  is  7.5-8s  after 
the  direct  P.  However,  at  all  three  stations,  we  have  not  observed  "stable"  Moho  Ps  on 
the  receiver  functions  profiles  (Fig.  5).  ERDO’s  receiver  functions  are  almost  dom¬ 
inated  by  the  converted  waves  from  shallow  structure  so  that  it  is  difficult  to  identify 
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the  Moho  Ps. 

We  found  that  the  receiver  function  is  an  efficient  way  to  study  the  lateral  varia¬ 
tion  of  crustal  structure,  provided  there  are  substantial  teleseismic  events  with  good 
back-azimuth  distribution  have  been  recorded.  Modeling  both  the  radial  and  tangential 
receiver  function  requires  ray-tracing  in  3D  medium.  Although  we  have  modeled  suc¬ 
cessfully  the  variation  pattern  of  radial  receiver  function  with  dipping  planar  interface, 
the  large  tangential  component  amplitudes  remain  hard  to  match.  More  sophisticated 
techniques  such  as  finite  difference  method  might  be  needed. 
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Table  1 

Back  Azimuthal/Distant  Ranges  and  Number  of  Events  of  Stacking  Suites 


Baz° 

TUNL 

A0 

number 

Baz° 

BUDO 

A0 

number 

Baz° 

ERDO 

A° 

number 

6±  1 

45±  1 

2 

5±  0 

46±  0 

1 

27±  1 

89±  5 

7 

27±  1 

90±  5 

5 

26±  1 

92±  5 

4 

45+  3 

54±  7 

10 

43±  4 

56±  8 

11 

41±  5 

59±10 

10 

58±  4 

40±  4 

36 

59±  6 

41±  4 

33 

57±  5 

42±  4 

20 

68±  4 

36±  1 

11 

78±  3 

36±  2 

14 

79+  5 

38±  2 

19 

75±  6 

39±  2 

9 

88±  5 

40±  3 

21 

89±  2 

42±  3 

7 

89±  4 

44±  4 

5 

102±  5 

48±  1 

5 

101±  4 

49±  1 

4 

no±  i 

84±  3 

9 

115±  4 

62±12 

14 

115±  4 

68±  7 

18 

112±  4 

71±10 

13 

127±  5 

46±  8 

20 

130±  6 

47±  7 

31 

128±  6 

48±  8 

27 

138±  5 

47±  6 

23 

145±  3 

50±  2 

3 

143±  3 

50±  2 

.  3 

153±  4 

48±  1 

4 

166±  1 

41±  1 

5 

161±  5 

43±  3 

5 

164±  0 

39±  0 

1 

182±  0 

52±  0 

2 

181±  0 

51±  0 

2 

180±  0 

50±  0 

2 

212±  0 

59±  0 

1 

211±  0 

58±  0 

1 

21G±  0 

57±  0 

1 

243±  0 

40±  0 

1 

255±  0 

59+7 

2 

255±  0 

51±  0 

1 

266±  0 

32±  0 

2 

272±  1 

37±  0 

2 

270±  4 

33±  2 

3 

272±  3 

33±  1 

2 

285±  0 

37±  2 

2 

290±  4 

46±  7 

6 

291±  1 

47±  5 

4 

289±  3 

38±  4 

4 

303±  0 

54±  0 

1 

303±  0 

53±  0 

1 

303±  0 

53±  0 

1 

346±  0 

54±  0 

1 

346±  0 

55±  0 

1 

Total 

205 

161 

102 
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Figure  Captions 


Figure  1.  Tectonic  sketch  map  of  Tibet  Plateau  (modified  from  Dewey  et  al.,  1988). 
Grayed  areas  are  below  3  km  elevation. 


Figure  2.  Synthetic  receiver  function  variation  with  back-azimuth  for  a  model  contain¬ 
ing  a  mid-crust  low  velocity  layer  with  its  upper  interface  dipping  20°  to  the  south  (see 
the  paper).  P-to-S  converted  phase  from  the  dipping  interface  is  indicated  by  dashed 
line. 


Figure  3.  (a)  Amplitude  variation  with  back-azimuth  of  P-to-S  converted  phase  from 
dipping  interfaces  at  different  dip  angles,  (b)  maximum  amplitude  variation  with  velo¬ 
city  contract  across  the  dipping  interface  at  different  dip  angles.  The  S-wave  velocity 
differences  are  from  0.3  to  0.8  km/s  correspond  to  S-wave  velocity  contrasts  of 
3.4/3. 1,  3.5/3. 1,  3.5/3.0,  3.6/3.0,  3.6/2.9,  and  3.6/2.8  across  the  interface.  Teleseismic 
ray  parameter  is  0.068  s/km. 


Figure  4.  Teleseismic  events  used  in  this  receiver  function  study,  their  epicentral  dis¬ 
tances  range  from  35°  to  95°  with  respect  to  station  BUDO.  The  area  in  figure  1  is 
shown  as  gray  box  in  the  figure. 


Figure  5.  Variation  of  receiver  functions  with  back-azimuth  in  station  (a)  TUNL,  (b) 
BUDO,  and  (c)  ERDO.  Arrows  point  to  the  back-azimuth  with  minimum  tangential 
energy.  Also  indicated  by  blackened  dots  are  mid-crust  P-to-S  converted  phases. 


Figure  6.  Measured  absolute  amplitudes  of  the  mid-crust  converted  phases  at  three  sta¬ 
tions.  Dashed  error  bar  means  only  one  event  is  available  in  this  stacking  suite  so 
average  standard  deviation  of  other  stacking  suites  is  used.  Solid  line  are  synthetic 
calculations  with  the  velocity  contrast  and  dip  angle  shown  in  the  figure. 


Figure  7.  Synthetic  receiver  functions  showing  the  effect  of  shallow  dipping  interface 
under  ERDO.  Dashed  line  indicats  the  Ps  phase  from  the  dipping  interface  (see  paper). 
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Figure  1.  Tectonic  sketch  map  of  Tibet  Plateau  (modified  from  Dewey  et  al.,  1988). 
Grayed  areas  are  below  3  km  elevation. 


Figure  2.  Synthetic  receiver  function  variation  with  back-azimuth  for  a  model  contain¬ 
ing  a  mid-crust  low  velocity  layer  with  its  upper  interface  dipping  20°  to  the  south  (see 
the  paper).  P-to-S  converted  phase  from  the  dipping  interface  is  indicated  by  dashed 
line. 


Figure  3.  (a)  Amplitude  variation  with  back-azimuth  of  P-to-S  converted  phase  from 
dipping  interfaces  at  different  dip  angles,  (b)  maximum  amplitude  variation  with  velo¬ 
city  contract  across  the  dipping  interface  at  different  dip  angles.  The  S-wave  velocity 
differences  are  from  0.3  to  0.8  km/s  correspond  to  S-wave  velocity  contrasts  of 
3.4/3.1,  3 .5/3.1,  3.5/3.0,  3.6/3.0,  3.6/2.9,  and  3.6/2.8  across  the  interface.  Teleseismic 
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Figure  6.  Measured  absolute  amplitudes  of  the  mid-cnist  converted  phases  at  three  sta¬ 
tions.  Dashed  error  bar  means  only  one  event  is  available  in  this  stacking  suite  so 
average  standard  deviation  of  other  stacking  suites  is  used.  Solid  line  are  synthetic 
calculations  with  the  velocity  contrast  and  dip  angle  shown  in  the  figure. 
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Figure  7.  Synthetic  receiver  functions  showing  the  effect  of  shallow  dipping  interface 
under  ERDO.  Dashed  line  indicats  the  Ps  phase  from  the  dipping  interface  (see  paper). 
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OBJECTIVE 

Broadband  records  of  regional  earthquakes  can  help  to  constrain  the  source  mechanisms 
and  crustal  structure.  The  objective  of  this  report  is  to  study  the  crustal  structure  and  source 
mechanisms  by  using  regional  earthquake  waveform  data  obtained  during  the  1991-1992  pas¬ 
sive  recording  experiment  on  the  Tibetan  Plateau  (Owens  et  al.,  1993).  The  data  also  provide 
an  opportunity  to  test  the  ability  of  determining  source  mechanisms  of  regional  events  by 
broadband  waveform  modeling  in  an  area  where  only  limited  number  of  stations  are  available 
and  the  structure  is  not  well  known. 


RESEARCH  ACCOMPLISHED 
Introduction 

The  Tibet  Plateau,  produced  by  collision  of  two  continental  plates,  India  and  Eurasia,  has 
been  a  challenge  in  Earth  science  for  a  long  time.  Several  models  of  collision  and  formation  of 
the  plateau  have  been  proposed.  But  none  of  them  can  give  a  satisfactory  explanation. 

Seismic  velocity  structure  of  the  crust  and  upper  mantle  can  provide  us  information  on 
the  present  state  deep  in  the  plateau,  which  is  crucial  to  understanding  the  process  of  plateau 
formation.  Due  to  the  difficulties  of  accessing  this  area,  previous  seismic  studies  had  to  rely  on 
the  data  recorded  at  stations  outside  the  plateau.  Most  studies  were  done  by  using  long-period 
surface  wave,  short-period  travel  time  analysis,  or  waveform  modeling  (S,  SS  etc),  see  Zhao  et 
al  (1991).  In  1991-1992,  a  passive  source  seismic  recording  experiment  was  conducted  on  the 
Tibet  Plateau  by  Sino-US  seismologists.  For  the  first  time,  11  digital  broadband  seismic  sta¬ 
tions  were  deployed  on  the  plateau  (Figure  1).  Large  amount  of  data  was  obtained  from  telese- 
ismic  and  regional  events. 

In  this  report,  we  will  use  the  data  from  several  regional  events  to  study  the  source 
mechanisms  and  the  crustal  structure.  Two  previous  velocity  models  are  examined  and  a  new 
model  is  formed  from  the  waveform  modeling.  We  will  also  test  the  ability  of  of  constraining 
the  source  mechanisms  of  regional  events  by  one  or  few  broadband  stations. 

Relocation  of  events  and  Pn,  Sn  velocities 

More  than  50  regional  events  with  reasonable  signal/noise  ratio  have  been  recorded  dur¬ 
ing  the  one  year  experiment.  We  selected  5  events  with  4  of  them  located  within  the  array 
(Figure  1).  Due  to  lack  of  reports  of  regional  seismic  stations  and  the  difference  of  Tibet  crust 
structure  from  standard  earth  model,  the  source  location  parameters  of  the  events  in  Tibet  give 
by  PDE  are  not  satisfactory.  Thus,  we  first  relocate  the  events  by  their  first  arrivals.  The  velo¬ 
city  model  for  relocation  is  a  two-layer  crust  model  simplified  from  previous  teleseismic 
receiver  function  study  (Zhu,  1993).  It  has  a  3  km  top  layer  with  P  velocity  4.5  km/s  and  62 
km  layer  with  P  velocity  of  6.2  km/s,  the  uppermost  mantle  velocity  is  8.1  km/s.  Since  most 
earthquakes  in  Tibet  occur  at  depths  between  5  to  20  km  (Molnar  and  Chen,  1983;  Zhao  and 
Helmberger  1991),  we  fix  the  source  depth  at  10  km.  The  relocation  results  are  listed  in  Table 
1.  It  is  shown  that  PDE  tends  to  give  later  origin  times  for  events  in  Tibet  plateau  because  the 
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crust  thickness  used  by  PDE  is  much  less  than  65  km  average  thickness  of  Tibet  crust. 


Table  1  Event  list  and  relocation  results 


PDE  n 

Relocation 

Event 

origin  time  lat.  long,  h  Mb  rms 
yr-mo-da  hr  mi  tO  0  0  km  s 

tO  lat.  long,  rms  AtO  Ax 
s  0  0  s  s  km 

222 

91-08-10  20  21  51.7  33.91  92.16  10  4.7  1.2 

52.1  33.92  92.28  0.9  0.4  11 

323 

91-11-19  01  04  18.0  32.48  93.59  33  4.9  3.0 

15.9  32.55  93.82  0.8  -2.1  23 

330 

91-11-26  21  15  59.9  34.07  94.25  33  4.3  2.5 

57.6  34.09  94.23  0.8  -2.3  3 

336 

91-12-02  19  45  36.6  32.09  94.69  *  4.4  4.0 

33.6  32.17  94.62  0.7  -3.0  11 

348 

91-12-14  08  20  23.8  33.98  88.84  33  5.1  5.4 

20.8  33.93  89.06  0.3  -3.0  20 

Event  222  was  recorded  by  all  the  11  ‘stations  in  epicentral  distance  range  of  80-800  km. 
In  addition,  it  occurred  in  the  middle  of  the  north-south  station  line,  forming  a  nature  seismic 
profile.  Figure  2  is  a  record  section  of  this  event.  By  identifying  the  Pn  and  Sn  arrivals  and 
using  least-square  fitting,  we  obtain  the  following  travel  time  relations: 

Pn  =  12.39+A/8.11  (sec) 

Sn  =  23.90+A/4.71  (sec) 

Our  8.11  km/s  and  4.71  km/s  for  P  and  S  velocities  in  uppermost  mantle  are  very  consistent 
with  previous  results  (Chen  and  Molnar,  1981)  The  average  crust  P  and  S  velocities  estimated 
from  intercept  time  are  3.46  km/s  and  6.24  km/s,  assuming  the  source  depth  of  10  km  and 
crust  thickness  of  65  km.  All  these  velocity  analyses,  combined  with  results  of  previous  work, 
provide  initial  models  for  our  next  step  involving  waveform  modeling. 

Waveform  modeling 

To  estimate  the  source  mechanisms,  we  used  a  procedure  developed  by  Zhao  and  Helm- 
berger  (1994).  It  uses  a  direct  grid  search  in  the  strike,  dip,  rake  and  source  depth  parameter 
space  for  the  minimum  LI  and  L2  norms  of  the  difference  between  wavefonns  and  synthetics. 
The  method  desensitizes  the  timing  between  principal  crustal  arrivals  by  fitting  portions  of  the 
waveforms  independently,  allowing  some  tune  shifts  between  the  observation  and  synthetics. 
We  find  that  introducing  one  more  parameter,  the  scalar  moment,  in  the  search  process  and 
using  the  same  moment  for  all  stations  and  components  increases  the  resolution. 

Two  previous  velocity  models  of  Tibet  plateau  have  been  examined.  Model  M45  was 
derived  from  pure  path  phase  velocity  dispersions  of  long-period  surface  wave  (Romanowicz, 
1982).  Another  model,  TIB,  was  from  waveform  modeling  of  regional  long-period  Love  wave 
(Zhao  et  al.,  1991).  Figure  3  shows  these  two  models.  Both  models  give  reasonably  good 
waveform  fits  of  the  Love  waves.  But  it  turns  out  that  TIB  yields  SmS  arrivals  that  are  too  fast 
relative  to  the  observed  SmS  (Figure  4).  M45’s  prediction  of  SmS  arrival  times  is  better  than 
TIB’s.  By  comparing  these  two  model,  we  conclude  that  TIB  has  mid-lower  crustal  velocities 
that  are  too  high.  This  conclusion  is  also  consistent  with  previous  receiver  function  study 
(Zhu,  1993). 

We  start  with  a  simple  initial  model:  a  low  velocity  top  layer  on  a  uniform  crustal  layer 
of  50  km  thick  with  a  transition  zone  above  moho.  By  adjusting  the  thickness  and  velocities  of 
each  layer,  we  try  to  model  the  observed  waveforms  of  both  Pnl  and  surface  wave.  Generally, 
the  Love  wave  provides  constraint  on  the  velocity  structure  of  upper  crust  and  Pnl,  SmS  pro¬ 
vide  constraint  on  the  lower  crust.  Our  final  model  is  given  in  Figure  3  labeled  T91:  the  top 
layer  thickness  is  3  km  with  a  S  velocity  of  2.50  km/s  and  a  P  velocity  of  4.45  km/s  followed 
by  a  50  km  main  crust  with  a  S  velocity  of  3.45  km/s  and  a  P  velocity  of  6.2  km/s.  Moho  is  at 
a  depth  of  68  km  with  a  15  km  transition  zone  above  it.  Figure  5  is  the  waveform  fit  by  this 
model.  For  all  5  events,  this  model  produces  good  fits  between  observations  and  synthetics, 
especially  the  SH  components.  The  obtained  source  mechanisms  are  listed  in  Table  2.  The  two 
northern  events  are  strike-slip  type.  The  other  3  are  normal  fault  type  (Figure  1). 
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Table  2  Source  mechanisms  from  waveform  modeling 


ent 

whole  array 
strike  dip  rake  ~  M0 

10zz  dyn  cm 

depth 

km 

error 

one  station  (AMDO) 
strike  dip  rake  --  M0  depth 

10z  dyn  cm  km 

error 

222 

340 

70 

170 

20.0 

10 

3.7 

340 

80  170 

21.0 

10 

3.2 

323 

190 

50 

260 

1.5 

10 

4.0 

170 

50  220 

1.8 

10 

3.1 

330 

30 

60 

210 

7.3 

10 

4.2 

30 

80  240 

8.8 

10 

3.4 

336 

170 

30 

220 

2.6 

10 

4.1 

210 

30  260 

2.3 

10 

3.3 

348 

0 

60 

240 

19.0 

5 

4.0 

190 

70  220 

20.0 

10 

3.3 

CONCLUSION  AND  RECOMMENDATIONS  - 

The  result  of  regional  earthquake  waveform  modeling  shows  that  the  crust  of  Tibet  pla¬ 
teau  is  roughly  vertically  uniform  and  has  a  very  low  average  velocity,  especially  the  shear 
velocity.  It  suggests  that  the  temperature  in  the  mid-lower  crust  is  high  which  might  be  caused 
by  the  concentration  of  continental  crust  materials  or  high  heat  flow  from  the  upper  mantle. 
Our  results  are  consistent  with  previous  surface  wave  and  teleseismic  receiver  function  studies. 
Two  types  of  source  mechanisms  have  been  obtained  for  5  regional  events:  strike-slip  with  one 
EW-striking  fault  plane  and  normal  fault  with  EW  T-axes.  The  source  mechanisms  and  then- 
indicated  stress  states  support  the  proposed  EW  extrusion  and  extension  of  the  crust  in  the  high 
plateau  under  the  indentation  of  India  subcontinent  (Molar  and  Layon-Caen,  1989). 

To  test  the  ability  of  constraining  source  mechanism  of  regional  events  by  one  or  few 
broadband  stations,  we  select  one  station  (AMDO)  and  only  use  its  recordings  to  search  for  the 
source  mechanisms.  The  solutions  are  also  listed  in  Table  2  with  solutions  from  the  whole 
array.  The  result  is  encouraging.  Under  the  circumstance  that  we  have  good  Pnl  and  surface 
wave  recordings  and  suitable  velocity  model,  the  source  mechanisms  can  be  obtained  by  one 
station,  as  found  in  California  (Dreger  and  Helmberger,  1993).  We  believe  such  kind  of  test  is 
very  helpful  for  further  broadband  waveform  modeling  in  area  as  Tibet  where  only  few  stations 
are  available.  For  example,  beginning  from  1992  CDSN  has  one  more  permanent  broadband 
station  operating  at  Lhasa.  More  regional  earthquake  data  can  be  accumulated  in  recent  future. 
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Tectonic  sketch  map  of  Tibet  plateau  (modified  from  Dewey  et  a!.,  1988),  gray  areas  are  below  3  km  elevation, 
Blacken  dots  are  5  events  with  their  source  mechanisms  from  waveform  modeling 
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Fig.  2  Recording  section  of  event  222,  dashed  lines  show  the  Pn,  Sn  and  SmS. 
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Fig.  5a  Waveform  modeling  of  Event  222 ,  upper  traces  are  observation  and  lower  traces  are 
synthetics  with  time  shifts  in  sec  (positive  number  means  synthetic  is  shifted  back).  Numbers 
below  station  names  are  distances  in  km.  Model  used  is  T91  with  source  at  depth  of  10  km. 
Source  mechanism  from  modeling  is  strike  340,  dip  70,  rake  170,  and  moment  2.0  x  1023  dyn. 
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Fig.  5b  Waveform  modeling  of  Event  323,  upper  traces  are  observation  and  lower  traces  are 
synthetics  with  time  shifts  in  sec  (positive  number  means  synthetic  is  shifted  back).  Numbers 
below  station  names  are  distances  in  km.  Model  used  is  T91  with  source  at  depth  of  10  km. 

Source  mechanism  from  modeling  is  strike  190,  dip  50,  rake  260,  and  moment  1.5  x  1022  dyn. 
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Fig.  5c  Waveform  modeling  of  Event  330,  upper  traces  are  observation  and  lower  traces  are 
synthetics  with  time  shifts  in  sec  (positive  number  means  synthetic  is  shifted  back).  Numbers 
below  station  names  are  distances  in  km.  Model  used  is  T91  with  source  at  depth  of  10  km. 
Source  mechanism  from  modeling  is  strike  30,  dip  60,  rake  210,  and  moment  7.3  x  1022  dyn.  cm 
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Fig.  Sd  Waveform  modeling  of  Event  336,  upper  traces  are  observation  and  lower  traces  are 
synthetics  with  time  shifts  in  sec  (positive  number  means  synthetic  is  shifted  back).  Numbers 
below  station  names  are  distances  in  km.  Model  used  is  T91  with  source  at  depth  of  10  km. 

Source  mechanism  from  modeling  is  strike  170,  dip  30,  rake  220,  and  moment  2.6  x  10^  dyn.  ^  40SCC 
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Fig.  5e  Waveform  modeling  of  Event  348,  upper  traces  are  observation  and  lower  traces  are 
synthetics  with  time  shifts  in  sec  (positive  number  means  synthetic  is  shifted  back).  Numbers 
below  station  names  are  distances  in  km.  Model  used  is  T91  with  source  at  depth  of  10  km. 

Source  mechanism  from  modeling  is  strike  10,  dip  50,  rake  250,  and  moment  2.0  x  1023  dyn.  cm 
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